############################################
# NOTE: The replication is conducted in the order that tables and figures are placed in paper.
############################################

############################################ ################################################
######################################## PACKAGES ##########################################
############################################ ################################################

if (!require("pacman")) {
  install.packages("pacman")
  library(pacman)
}

pacman::p_load(tidyverse, pwr, ggplot2, janitor, lubridate, readstata13, countrycode,
               estimatr, rdrobust, entropy, patchwork, abind, AER, aod, arsenal,
               bookdown, broom, caret, cobalt, corrplot, data.table, devtools, did2s,
               digest, dplyr, dygraphs, evtree, fabricatr, fBasics, foreach,
               foreign, future, geosphere, GGally, gganimate, ggExtra, ggmap, ggraph,
               ggstatsplot, gifski, glue, grf, grid, gridExtra, haven, hrbrthemes, jtools,
               kableExtra, cowplot, knitr, leaflet, lmtest, magrittr, MatchIt,
               memisc, modelsummary, multiwayvcov, openxlsx, sjlabelled, sjmisc, sp, stringr,
               texreg, tibble, tidyquant, tidyr, tidytext, xtable, rio, sjPlot, margins,
               sandwich, ggeffects, WeightIt)

############################################ ################################################
########################### SCRIPT, WORKING DIRECTORY & DATA ##########################################
############################################ ################################################

# function to create folders if they don't exist
if (!dir.exists("data")) {
  dir.create("data")
}

if (!dir.exists("out")) {
  dir.create("out")
}

if (!dir.exists("scripts")) {
  dir.create("scripts")
}

# clean environment
rm(list = ls())

# set your own working directory here
setwd("~/Dropbox/Russia Ukraine EU/Replication") 

# main dataset
dat = read_rds("data/dat.rds") 

# script log
sink("scripts/analysis_script_log.txt", append = TRUE) 

############################################ ################################################
######################################## MAIN TEXT ##########################################
############################################ ################################################

############################################ Figure 1 ############################################

dat$Treatment=0
dat$Treatment[dat$Date > as.Date("2022-02-23")]=1

dep_vars <- c("eu_support_st")
var_desc <- c("No better future outside the EU")
models <- list()

model_basic <- lm_robust(eu_support_st ~ Treatment, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_extended <- lm_robust(eu_support_st ~ Treatment + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_full <- lm_robust(eu_support_st ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_balanced <- lm_robust(eu_support_st ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)

models[["Basic model"]] <- model_basic
models[["Extended model"]] <- model_extended
models[["Full model"]] <- model_full
models[["Balanced model"]] <- model_balanced

p1 <- modelplot(models, coef_omit = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)") +
  labs(y = 'Treatment',
       title = "Support for the EU", subtitle = element_text(face="italic","EU sample")) + geom_vline(xintercept = 0, color = 'grey')+
  coord_flip() + 
  theme(legend.position = "bottom", legend.justification="left") +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        text = element_text(size = 12))
p1
ggsave(plot = p1,
       filename = "out/result_eu_st.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Figure 2 ############################################

days <- c(7, 14, 21, 30)
dat_list <- lapply(days, function(d) {
  dat %>%
    filter(day <= d)
})

dat7 <- dat_list[[1]]
dat14 <- dat_list[[2]]
dat21 <- dat_list[[3]]
dat30 <- dat_list[[4]]

model_rc7<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat7, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_rc14<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat14, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_rc21<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat21, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_rc30<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat30, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)

texreg(list(model_rc7, model_rc14, model_rc21, model_rc30), 
       booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE,
       caption = "Over-time Effects", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-overtime",
       float.pos = "h",
       custom.gof.rows = list("Region FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("Yes", "Yes", "Yes", "Yes"),
                              "Covariates" = c("Yes", "Yes", "Yes", "Yes"),
                              "Weight" = c("Yes", "Yes", "Yes", "Yes")))

models <- list(
  "1 week" = model_rc7,
  "2 weeks" = model_rc14,
  "3 weeks" = model_rc21,
  "4 weeks (Full)" = model_rc30)

p2 <- modelplot(models, coef_omit = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)")+
  labs(x = 'Coefficient estimates and 95% confidence intervals', 
       y = 'Term names',
       title = "Support for the EU", 
       subtitle = "Treatment effects over time") +
  guides(color = guide_legend(reverse = TRUE))  + geom_vline(xintercept = 0, color = 'grey') + coord_flip()
p2


ggsave(plot = p2,
       filename = "out/overtime_plot.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Figure 3 ############################################

# NOTE: Please note that figure 3 is replicated in the separate R script called forest_replication.R

############################################ ################################################
######################################## APPENDIX ###########################################
############################################ ################################################

############################################ Table 1 ############################################

treat_eu <- tabyl(dat, Treatment)
treat_country <- tabyl(dat, country_name, Treatment) %>% 
  mutate(n = `0` + `1`, treat_share = round(`1`/n, 2))
names(treat_country) <- c("Country", "Control", "Treated", "Total", "Treat Share")
knitr::kable(treat_country)
treat_country2 <- tabyl(dat, country_name, Treatment) %>% 
  mutate(n = `0` + `1`)
kable(treat_country2, booktabs = TRUE)
write.table(treat_country2, file = "out/Table1.txt", sep = "\t", quote = FALSE, row.names = FALSE)

############################################ Table 2 ############################################

desc <- dat %>% dplyr::select(age, female, did_not_attend_school, primary_school,
                              sec_grades_6_7, sec_grades_1_5, national_diploma,
                              higher_national_diploma, bachelors_degree, masters_degree,
                              doctoral_degree, education_level_10, education_level_11,
                              education_level_12, education_level_13, education_level_14,
                              education_level_15, education_level_16, education_level_17,
                              education_level_18, other_education, self_employed, managers,
                              other_white_collar, manual_workers, house_persons, unemployed,
                              retired, students, married_remarried, married_single_partner,
                              single, divorced_separated, widow, other_marital_status,
                              rural, small_town, large_town, internetuse, non_eu, ideology_num, Treatment)

balance_test_summary <- summary(tableby(Treatment ~ ., data = desc), title = "Balance Test by Treatment", text = TRUE)
capture.output(balance_test_summary, file = "out/Table2 and 3.txt")

############################################ Table 3 ############################################

# NOTE: The code above produces both tables but they are divided into two in Overleaf due to space constraints.

############################################ Table 4 ############################################

desc = dat %>% dplyr::select(eu_support, Treatment)
summary(tableby(Treatment ~ ., data = desc), title = "Dependent Variable by Treatment", text="latex")
dependent_variable_summary <- summary(tableby(Treatment ~ ., data = desc), title = "Dependent Variable by Treatment", text = TRUE)
capture.output(dependent_variable_summary, file = "out/Table4.txt")

############################################ Figure 4 ############################################

desc = dat %>% dplyr::select(Treatment, country_name, eu_support)  %>%
  group_by(country_name, Treatment) %>%
  summarise(mean_eu_support = round(mean(eu_support, na.rm = TRUE), 5))

desc <- desc %>%
  mutate(change_in_means = mean_eu_support[Treatment==1] - mean_eu_support[Treatment==0])

mean_row <- dat %>% 
  group_by(Treatment) %>%
  summarise(mean_eu_support = round(mean(eu_support, na.rm = TRUE), 5))

mean_row$country_name = "EU"

mean_row <- mean_row %>%
  mutate(change_in_means = mean_eu_support[Treatment==1] - mean_eu_support[Treatment==0])

descc <- bind_rows(desc, mean_row)

descc$Treatment[descc$Treatment==0]="No"
descc$Treatment[descc$Treatment==1]="Yes"

p4 = ggplot(descc, aes(x=mean_eu_support, y=reorder(country_name, mean_eu_support), fill=Treatment)) +
  geom_bar(stat="identity", position=position_dodge(width=0.5)) +
  theme_minimal() +
  labs(x="Mean EU Support", y="Country") +
  scale_fill_manual(values=c("black", "peachpuff2"), breaks=c("Yes", "No"))+
  theme_classic() +
  theme(plot.title = element_text(size=25),
        plot.subtitle = element_text(size=14))
p4

ggsave(plot = p4,
       filename = "out/mean_support_absolute.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Table 5 ############################################

dat7_summary <- dat7 %>% 
  group_by(Treatment) %>% 
  summarize(Count = n())

dat14_summary <- dat14 %>% 
  group_by(Treatment) %>% 
  summarize(Count = n())

dat21_summary <- dat21 %>% 
  group_by(Treatment) %>% 
  summarize(Count = n())

dat30_summary <- dat30 %>% 
  group_by(Treatment) %>% 
  summarize(Count = n())

print("Summary for day 7:")
print(dat7_summary)

print("Summary for day 14:")
print(dat14_summary)

print("Summary for day 21:")
print(dat21_summary)

print("Summary for day 30:")
print(dat30_summary)

combined_summary <- bind_rows(dat7_summary, dat14_summary, dat21_summary, dat30_summary)
write.table(combined_summary, file = "out/Table5.txt", sep = "\t", quote = FALSE, row.names = FALSE)

############################################ Figure 5 ############################################

p5 <- ggplot(dat, aes(x=Date)) + 
  geom_bar() +
  geom_vline(xintercept = as.Date("2022-02-24",format="%Y-%m-%d"),
             color =  "red",
             size = 0.5,
             linetype = "dashed") +
  labs(title = "Interviews per day before and after treatment",
       subtitle = "Red line indicates Russia's invasion of Ukraine on February 24, 2022",
       y = "Participants",
       x = "Day") +
  theme_classic() +
  theme(plot.title = element_text(size=25),
        plot.subtitle = element_text(size=14))
p5

ggsave(plot = p5,
       filename = "out/timing_EU.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Figure 6 ############################################

p6 <- ggplot(dat, aes(x=Date)) + 
  geom_bar() +
  geom_vline(xintercept = as.Date("2022-02-24",format="%Y-%m-%d"),
             color =  "red",
             size = 0.5,
             linetype = "dashed") +
  labs(title = "Interviews per day before and after treatment",
       subtitle = "Red line indicates Russia's invasion of Ukraine on February 24, 2022",
       y = "Participants",
       x = "Day") +
  theme_classic() +
  theme(plot.title = element_text(size=25),
        plot.subtitle = element_text(size=14)) +
  facet_wrap(~ country_name)
p6

ggsave(plot = p6,
       filename = "out/timing_countries.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Figure 7 ############################################
# Load data
google <- read_delim("data/multiTimeline.csv", delim = ";", escape_double = FALSE, trim_ws = TRUE)
google$Date <- as.Date(google$Day, "%d.%m.%Y")

p7 <- ggplot(google, aes(x = Date)) + 
  geom_line(aes(y = `ukraine: (Worldwide)`, color = 'Ukraine'), size = 0.5) + 
  geom_line(aes(y = `russia: (Worldwide)`, color = 'Russia'), size = 0.5) + 
  geom_vline(xintercept = as.numeric(as.Date("2022-02-24")), color = "grey", linetype = "dashed", size = 0.4) + 
  theme_minimal() + 
  labs(title = "",
       subtitle = "",
       y = "Search",
       x = "Date",
       color = "Group") + 
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5),
        axis.title = element_text(size = 10),
        plot.title = element_text(size = 12),
        plot.subtitle = element_text(size = 10),
        legend.title = element_text(size = 10),
        legend.position = "right") +
  scale_color_manual(values = c('Ukraine' = '#6CA6CD', 'Russia' = '#9370DB')) + 
  scale_x_date(date_labels = "%d-%b", date_breaks = "3 days", limits = as.Date(c("2022-02-21", "2022-03-23")))

p7

ggsave(plot = p7,
       filename = "out/google.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Figure 8 ############################################

n1 <- 1887   
n2 <- 22611
alpha_05 <- 0.05
alpha_10 <- 0.1
power_values <- seq(0.6, 1, by = 0.01)

effect_sizes_05 <- numeric(length(power_values))
effect_sizes_10 <- numeric(length(power_values))

for (i in seq_along(power_values)) {
  effect_sizes_05[i] <- pwr.t2n.test(n1 = n1, n2 = n2, sig.level = alpha_05, power = power_values[i])$d
  effect_sizes_10[i] <- pwr.t2n.test(n1 = n1, n2 = n2, sig.level = alpha_10, power = power_values[i])$d
}

data_for_plot <- data.frame(Power = rep(power_values, 2),
                            EffectSize = c(effect_sizes_05, effect_sizes_10),
                            Alpha = factor(rep(c("0.05", "0.1"), each = length(power_values))))

p8 = ggplot(data_for_plot, aes(x = Power, y = EffectSize, color = Alpha)) +
  geom_line() +
  theme_minimal() +
  scale_x_continuous(limits = c(0.6, 0.9)) + 
  scale_y_continuous(limits = c(0.045, 0.08)) + 
  labs(title = "Effect Size by Power",
       x = "Power",
       y = "Expected Effect Size (Cohen's d)",
       color = "alpha") +
  theme(legend.title = element_text(size = 8))

print(p8)

ggsave(plot = p8,
       filename = "out/power_ukraine.eps",
       width = 15,
       height = 12,
       units = "in")

############################################ Table 6 ############################################

model_names <- c("Basic model", "Extended model", "Full model", "Balanced model")

dat$cluster = as.factor(paste(dat$Date, dat$country_code, sep = "_"))

model_basic <- lm_robust(eu_support ~ Treatment, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_extended <- lm_robust(eu_support ~ Treatment + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_full <- lm_robust(eu_support ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_balanced <- lm_robust(eu_support ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)

texreg(list(model_basic, model_extended, model_full, model_balanced), 
       booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
       custom.model.names = model_names,
       caption = paste("Treatment and EU Support - EU sample"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = paste("tab:estimates-eulevel"),
       float.pos = "h",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

screenreg(list(model_basic, model_extended, model_full, model_balanced), 
          file = "out/Table6.txt",
          booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
          custom.model.names = model_names,
          caption = paste("Treatment and EU Support - EU sample"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
          omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", 
          custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                                 "Time" = c("-", "Yes", "Yes", "Yes"),
                                 "Covariates" = c("-", "-", "Yes", "Yes"),
                                 "Weight" = c("-", "-", "-", "Yes")))


############################################ Table 7 ############################################

model_names <- c("Basic model", "Extended model", "Full model", "Balanced model")

model_basic <- lm_robust(eu_support_st ~ Treatment, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_extended <- lm_robust(eu_support_st ~ Treatment + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_full <- lm_robust(eu_support_st ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster)
model_balanced <- lm_robust(eu_support_st ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)

texreg(list(model_basic, model_extended, model_full, model_balanced), 
       booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
       custom.model.names = model_names,
       caption = paste("Treatment and EU Support - EU sample"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = paste("tab:estimates-eulevel-st"),
       float.pos = "h",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

screenreg(list(model_basic, model_extended, model_full, model_balanced), 
          file = "out/Table7.txt",
          booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
          custom.model.names = model_names,
          caption = paste("Treatment and EU Support - EU sample"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
          omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", 
          custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                                 "Time" = c("-", "Yes", "Yes", "Yes"),
                                 "Covariates" = c("-", "-", "Yes", "Yes"),
                                 "Weight" = c("-", "-", "-", "Yes")))

############################################ Table 8 ############################################

model_rc7<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat7, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_rc14<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat14, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_rc21<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat21, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_rc30<-lm_robust(eu_support_st~Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day, data=dat30, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)

texreg(list(model_rc7, model_rc14, model_rc21, model_rc30), 
       booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE,
       caption = "Over-time Effects", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-overtime",
       float.pos = "h",
       custom.gof.rows = list("Region FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("Yes", "Yes", "Yes", "Yes"),
                              "Covariates" = c("Yes", "Yes", "Yes", "Yes"),
                              "Weight" = c("Yes", "Yes", "Yes", "Yes")))

screenreg(list(model_rc7, model_rc14, model_rc21, model_rc30), 
       booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE,
       file = "out/Table8.txt",
       caption = "Over-time Effects", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-overtime",
       float.pos = "h",
       custom.gof.rows = list("Region FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("Yes", "Yes", "Yes", "Yes"),
                              "Covariates" = c("Yes", "Yes", "Yes", "Yes"),
                              "Weight" = c("Yes", "Yes", "Yes", "Yes")))

############################################ Table 9 ############################################
dat$cluster = as.factor(paste(dat$Date, dat$region, sep = "_"))

dat$Borders_Russia=0
dat$Borders_Russia[dat$isocntry=="LV"]=1
dat$Borders_Russia[dat$isocntry=="EE"]=1
dat$Borders_Russia[dat$isocntry=="LT"]=1
dat$Borders_Russia[dat$isocntry=="PL"]=1
dat$Borders_Russia[dat$isocntry=="FI"]=1
dat$Borders_Russia <- haven::labelled(dat$Borders_Russia, label = "Bordering Russia")

models <- list()
tables_list <- list()
dep_vars <- c("eu_support_st")
var_desc <- c("No better future outside the EU")
model_names <- c("Basic", "Extended", "Full", "Balanced")

for (v in dep_vars) {
  
  model_basic <- lm_robust(as.formula(paste(v, "~ Treatment*Borders_Russia")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster)
  model_extended <- lm_robust(as.formula(paste(v, "~ Treatment*Borders_Russia + day")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster)
  model_full <- lm_robust(as.formula(paste(v, "~ Treatment*Borders_Russia + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster)
  model_balanced <- lm_robust(as.formula(paste(v, "~ Treatment*Borders_Russia + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster, weights = webal)
  
  models[[v]] <- list(model_basic, model_extended, model_full, model_balanced)
  
  print(texreg(list(model_basic, model_extended, model_full, model_balanced), 
               booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
               custom.model.names = model_names,
               caption = paste("Treatment and EU Support - HTE by bordering Russia"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
               omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = paste("tab:estimates-", v ,"-borderlevel"),
               float.pos = "h",
               custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                                      "Time" = c("-", "Yes", "Yes", "Yes"),
                                      "Covariates" = c("-", "-", "Yes", "Yes"),
                                      "Weight" = c("-", "-", "-", "Yes"))))
}

screenreg(list(model_basic, model_extended, model_full, model_balanced), 
             booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE,
             custom.model.names = model_names,
          file = "out/Table9.txt",
             caption = paste("Treatment and EU Support - HTE by bordering Russia"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
             omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = paste("tab:estimates-", v ,"-borderlevel"),
             float.pos = "h",
             custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                                    "Time" = c("-", "Yes", "Yes", "Yes"),
                                    "Covariates" = c("-", "-", "Yes", "Yes"),
                                    "Weight" = c("-", "-", "-", "Yes")))

############################################ Table 10 ###########################################

bordering_countries <- c("LV", "EE", "LT", "PL", "FI", "SK", "HU", "RO")
dat$Borders_Russia2 <- ifelse(dat$isocntry %in% bordering_countries, 1, 0)
dat$Borders_Russia2 <- haven::labelled(dat$Borders_Russia2, label = "Bordering Russia or Ukraine")

eastern_europe_countries <- c("BG", "CY", "CZ", "EE", "GR", "HR", "HU", "LT", "LV", "PL", "RO", "SI", "SK")
eastern_europe_countries_and_finland <- c("BG", "CY", "CZ", "EE", "GR", "HR", "HU", "LT", "LV", "MT", "PL", "RO", "SI", "SK", 'FI')

dat$eastern_europe <- ifelse(dat$isocntry %in% eastern_europe_countries, 1, 0)
dat$eastern_europe_finland <- ifelse(dat$isocntry %in% eastern_europe_countries_and_finland, 1, 0)

W2<- weightit(Treatment ~ age + female + did_not_attend_school + primary_school +
                sec_grades_6_7 + sec_grades_1_5 + national_diploma +
                higher_national_diploma + bachelors_degree + masters_degree +
                doctoral_degree + education_level_10 + education_level_11 +
                education_level_12 + education_level_13 + education_level_14 +
                education_level_15 + education_level_16 + education_level_17 +
                education_level_18 + other_education + self_employed + managers +
                other_white_collar + manual_workers + house_persons + unemployed +
                retired + students + married_remarried + married_single_partner +
                single + divorced_separated + widow + other_marital_status +
                rural + small_town + large_town + internetuse + non_eu + ideology_num + eastern_europe, data = dat, tmethod = "ebal", estimand = "ATE", by = "country_name")
webal_c <- W2$weights
datt <- cbind(dat, webal_c)

models <- list()
tables_list <- list()
dep_vars <- c("eu_support_st")
var_desc <- c("No better future outside the EU")
model_names <- c("Balanced", 'Balanced', "Balanced")

for (v in dep_vars) {
  
  model_balanced0 <- lm_robust(as.formula(paste(v, "~ Treatment*Borders_Russia2 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster, weights = webal)
  model_balanced1 <- lm_robust(as.formula(paste(v, "~ Treatment*eastern_europe + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster, weights = webal)
  model_balanced2 <- lm_robust(as.formula(paste(v, "~ Treatment*eastern_europe_finland + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(region), cluster = cluster, weights = webal)
  
  models[[v]] <- list(model_balanced0,model_balanced1, model_balanced2)
  
  print(texreg(list(model_balanced0,model_balanced1, model_balanced2), 
               booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
               custom.model.names = model_names,
               caption = paste("Treatment and EU Support - HTE by Eastern Europe"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
               omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)",
               label = paste("tab:estimates-", v ,"-easterneurope"),
               float.pos = "h",
               custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes"),
                                      "Time" = c("Yes", "Yes", "Yes"),
                                      "Covariates" = c("Yes", "Yes", "Yes"),
                                      "Weight" = c("Yes", "Yes", "Yes"))))
}

screenreg(list(model_balanced0,model_balanced1, model_balanced2), 
             booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
             custom.model.names = model_names,
          file = "out/Table10.txt",
             caption = paste("Treatment and EU Support - HTE by Eastern Europe"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
             omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)",
             label = paste("tab:estimates-", v ,"-easterneurope"),
             float.pos = "h",
             custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes"),
                                    "Time" = c("Yes", "Yes", "Yes"),
                                    "Covariates" = c("Yes", "Yes", "Yes"),
                                    "Weight" = c("Yes", "Yes", "Yes")))

############################################ Table 11 ###########################################

dat$cluster = as.factor(paste(dat$Date, dat$country, sep = "_"))

control_group <- dat[dat$Treatment == 0, ]
if (all(is.na(control_group))) {
  stop("Control group contains only NAs")
}

standardize <- function(var) {
  mean_var <- mean(control_group[[var]], na.rm = TRUE)
  sd_var <- sd(control_group[[var]], na.rm = TRUE)
  
  if (is.na(mean_var) || is.na(sd_var) || sd_var == 0) {
    warning(paste("Mean or SD is NA or SD is 0 for variable", var))
    return(rep(NA, nrow(dat)))
  }
  
  return((dat[[var]] - mean_var) / sd_var)
}

dat <- dat %>%
  mutate(national_democracy = case_when(
    sd18a == 1 ~ 4,
    sd18a == 2 ~ 3,
    sd18a == 3 ~ 2,
    sd18a == 4 ~ 1,
    TRUE ~ NA_real_
  )) %>%
  mutate(national_democracy = (national_democracy - 1) / 3)

dat$national_democracy_st <- standardize('national_democracy')

dat <- dat %>%
  mutate(national_right_direction = case_when(
    d73_1 == 1 ~ 1,
    d73_1 == 2 ~ 0,
    d73_1 == 3 ~ NA_real_,
    TRUE ~ NA_real_
  ))
dat$national_right_direction_st <- standardize('national_right_direction')

models <- list()
dep_vars <- c("national_democracy_st", "national_right_direction_st")
model_names <- c("Balanced")

model_balanced_national_democracy <- lm_robust(national_democracy_st ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
model_balanced_national_right_direction <- lm_robust(national_right_direction_st ~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)

models[[1]] <- model_balanced_national_democracy
models[[2]] <- model_balanced_national_right_direction

print(texreg(models, 
             booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
             caption = "Treatment and National Variables - EU sample", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
             omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", 
             label = "tab:treatment-national-variables-eu-sample", float.pos = "h",
             custom.gof.rows = list("Country FE" = c("Yes", 'Yes'),
                                    "Time" = c("Yes", 'Yes'),
                                    "Covariates" = c("Yes", 'Yes'),
                                    "Weight" = c("Yes", 'Yes'))))

screenreg(models, 
             booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
             caption = "Treatment and National Variables - EU sample", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
             omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", 
             label = "tab:treatment-national-variables-eu-sample", float.pos = "h",
       file = "out/Table11.txt",
             custom.gof.rows = list("Country FE" = c("Yes", 'Yes'),
                                    "Time" = c("Yes", 'Yes'),
                                    "Covariates" = c("Yes", 'Yes'),
                                    "Weight" = c("Yes", 'Yes')))

############################################ Figure 9 ###########################################
############################################ Figure 9 ###########################################

W <- weightit(Treatment ~ + female + did_not_attend_school + primary_school +
                sec_grades_6_7 + sec_grades_1_5 + national_diploma +
                higher_national_diploma + bachelors_degree + masters_degree +
                doctoral_degree + education_level_10 + education_level_11 +
                education_level_12 + education_level_13 + education_level_14 +
                education_level_15 + education_level_16 + education_level_17 +
                education_level_18 + other_education + self_employed + managers +
                other_white_collar + manual_workers + house_persons + unemployed +
                retired + students + married_remarried + married_single_partner +
                single + divorced_separated + widow + other_marital_status +
                rural + small_town + large_town + internetuse + non_eu + ideology_num, data = dat, tmethod = "ebal", estimand = "ATE", by = "country_name")
webal_c <- W$weights
datt <- cbind(dat, webal_c)

countries <- datt %>% pull(country_name) %>% unique()
get.res <- function(x){
  
  est <- x$coefficients["Treatment"]
  se <- x$std.error["Treatment"]
  ci_low <- x$conf.low["Treatment"]
  ci_high <- x$conf.high["Treatment"]
  
  res <- cbind(est, se, ci_low, ci_high)
  
  return(res)
  
}

model_name <- c("Balanced")
model_formula <- c("Treatment + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + day")
dep_var <- c("eu_support_st")
var_desc <- c("No better future outside the EU")
res_df <- data.frame()

for (v in 1:length(dep_var)) {
  
  for (i in 1:length(countries)){
    
    tmp_df <- datt %>% filter(country_name == countries[i])
    
    for (m in 1:length(model_name)){
      
      reg_form <- as.formula(paste(dep_var[v], "~", model_formula[m]))
      
      if(m == 1){res_reg <- lm_robust(reg_form,
                                      fixed_effects = ~ factor(region), cluster = Date,
                                      data = tmp_df, `se_type` = "stata", weights = webal_c)}
      tmp_res_df <- get.res(res_reg) %>%
        as.data.frame() %>%
        mutate(reg_type = model_name[m],
               country_name = countries[i],
               depvar = dep_var[v])
      res_df <- bind_rows(res_df, tmp_res_df)
      
    }
    
  }
  
}

custom_order <- rev(c("Greece", "Poland", "Slovakia", "Lithuania", "Portugal", "Cyprus", 
                      "Netherlands", "Luxembourg", "Croatia", "Hungary", "Romania", 
                      "Germany", "Slovenia", "France", "Sweden", "Malta", "Finland", 
                      "Belgium", "Denmark", "Italy", "Czechia", "Spain", "Estonia", 
                      "Ireland", "Latvia", "Austria", "Bulgaria"))

plot_list <- list()

for (v in 1:length(dep_var)) {
  
  plot_list[[v]] <- ggplot(data = filter(res_df, depvar == dep_var[v])) +
    geom_pointrange(mapping = aes(x = factor(country_name, levels = custom_order),
                                  y = est,
                                  ymin = ci_low,
                                  ymax = ci_high,
                                  colour = reg_type),
                    position = position_dodge(width = 0.8)) +
    coord_flip() +
    geom_hline(yintercept = 0, color = "grey", linetype = "dashed", size = 0.5) +
    theme_classic() +
    labs(title = "Effect of Russian Invasion of Ukraine on Feb 24, 2022",
         subtitle = "Dependent variable: EU Support",
         y = "Linear regression estimation with 95%-CI", 
         x = "Countries",
         colour = "Model specification") +
    theme(legend.position = "bottom") +
    scale_color_discrete(breaks = model_name)
  
  ggsave(plot = plot_list[[v]], 
         filename = paste0("out/", dep_var[v], "_plot.eps"),
         width = 12,
         height = 12,
         units = "in")
}
plot_list

############################################ Table 12 ###########################################

dat$eu_support=NA
dat$eu_support[dat$d19_2==1]= 1
dat$eu_support[dat$d19_2==2]= 2
dat$eu_support[dat$d19_2==3]= 3
dat$eu_support[dat$d19_2==4]= 4

dat$Treatment <- factor(dat$Treatment, levels = c(0, 1), labels = c("Control", "Treatment"))

dat$Ideology <- factor(ifelse(dat$ideology_num < 5, "Left", ifelse(dat$ideology_num > 5, "Right", "Centre")), levels = c("Centre", "Left", "Right"))
model <- lm_robust(eu_support ~ Treatment * Ideology + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, 
                   data = dat, 
                   se_type = "stata", 
                   fixed_effects = ~ factor(country), 
                   cluster = cluster, 
                   weights = webal)

print(texreg(list(model), 
             booktabs = TRUE, 
             dcolumn = TRUE, 
             include.ci = FALSE, 
             use.packages = FALSE, 
             custom.model.names = "Balanced",
             caption = "Treatment and EU Support - HTE by Ideological Categories", 
             stars = c(0.001, 0.01, 0.05), 
             include.rmse = FALSE, 
             omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", 
             label = "tab:estimates-ideology",
             float.pos = "h",
             custom.gof.rows = list("Country FE" = c("Yes"), "Time" = c("Yes"), "Covariates" = c("Yes"), "Weight" = c("Yes"))))

screenreg(list(model), 
             booktabs = TRUE, 
             dcolumn = TRUE, 
             include.ci = FALSE, 
             use.packages = FALSE, 
             custom.model.names = "Balanced",
             caption = "Treatment and EU Support - HTE by Ideological Categories",
          file = "out/Table12.txt",
             stars = c(0.001, 0.01, 0.05), 
             include.rmse = FALSE, 
             omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", 
             label = "tab:estimates-ideology",
             float.pos = "h",
             custom.gof.rows = list("Country FE" = c("Yes"), "Time" = c("Yes"), "Covariates" = c("Yes"), "Weight" = c("Yes")))

########################################### Figure 10 ###########################################

model <- lm(eu_support ~ Ideology*Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + factor(country), 
            data = dat)
clustered_se <- vcovHC(model, type = "HC1", cluster = ~dat$cluster)
coeftest(model, vcov = clustered_se)
ideology = plot_model(model, type = "pred", terms = c("Treatment", "Ideology"), 
                      title = "Marginal effects by ideology") +
  theme_minimal() +
  ylab("EU support")
ideology
ggsave("out/predicted_ideology_eu_support.eps", plot = ideology, device = "eps")

########################################### Table 13 ###########################################

dat$age_group <- cut(dat$age, 
                     breaks = c(0, 30, 40, 50, 100),  # Specify the age group breakpoints
                     labels = c("-30", "31-40", "41-50", "51+"),  # Labels for the age groups
                     include.lowest = TRUE)  # Include the lower bound in the first group
models <- list()
tables_list <- list()
dep_vars <- c("eu_support")
var_desc <- c("No better future outside the EU")
model_names <- c("Balanced")

for (v in dep_vars) {
  
  model_age <- lm_robust(as.formula(paste(v, "~ Treatment*age_group + day + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(country), cluster = cluster, weights = webal)
    models[[v]] <- list(model_age)
  
  print(texreg(list(model_age), 
               booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
               custom.model.names = model_names,
               caption = paste("Treatment and EU Support - HTE by age categories"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
               omit.coef = "day|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = paste("tab:estimates-", v ,"-age"),
               float.pos = "h",
               custom.gof.rows = list("Country FE" = c("Yes"),
                                      "Time" = c("Yes"),
                                      "Covariates" = c("Yes"),"Weight" = c("Yes"))))
}

screenreg(list(model_age), 
             booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE, use.packages = FALSE, 
             custom.model.names = model_names,
          file = "out/Table13.txt",
             caption = paste("Treatment and EU Support - HTE by age categories"), stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
             omit.coef = "day|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = paste("tab:estimates-", v ,"-age"),
             float.pos = "h",
             custom.gof.rows = list("Country FE" = c("Yes"),
                                    "Time" = c("Yes"),
                                    "Covariates" = c("Yes"),"Weight" = c("Yes")))

########################################### Figure 11 ###########################################

model <- lm(eu_support ~ Treatment * age_group + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu + factor(country), 
            data = dat)
clustered_se <- vcovHC(model, type = "HC1", cluster = ~dat$cluster)
coeftest(model, vcov = clustered_se)
age_plot <- plot_model(model, type = "int", terms = c("Treatment", "age_group"), 
                       title = "Marginal effects by age group") +
  theme_minimal() +
  ylab("EU Support")
age_plot
ggsave("out/predicted_age_group_eu_support.eps", plot = age_plot, device = "eps")

########################################### Figure 12 ###########################################
dat$cluster = as.factor(paste(dat$Date, dat$country, sep = "_"))

# NOTE: We need to code Treatment again and drop labels
dat$Treatment=0
dat$Treatment[dat$Date > as.Date("2022-02-23")]=1

vars_to_standardize <- c('eu_positive', 'eu_future_optimistic', 'eu_democracy', 'eu_right_direction', 'eu_political_matters', 'eu_future_outside', 'eu_more_decision', 'eu_my_voice', 'eu_country_voice', 'national_democracy', 'national_my_voice', 'national_political_matters', 'local_political_matters', 'national_right_direction')

for (var in vars_to_standardize) {
  if (all(is.na(dat[[var]]))) {
    warning(paste(var, "contains only NAs"))
  } else {
    dat[[paste0(var, "_st")]] <- standardize(var)
  }
}

dep_vars <- c("eu_positive_st",
              "eu_democracy_st",
              "eu_right_direction_st",
              "eu_more_decision_st",
              "eu_political_matters_st",
              "eu_future_optimistic_st",
              "eu_future_outside_st",
              "eu_my_voice_st",
              "eu_country_voice_st")

var_desc <- c("Positive image of EU",
              "Democratic satisfaction with EU",
              "EU in right direction",
              "More decision at the EU level",
              "Political discussion of EU",
              "Optimistic about future of EU",
              "No better future outside EU",
              "My voice counts in EU",
              "My country counts in EU")

models <- list()

for (i in seq_along(dep_vars)) {
  v <- dep_vars[i]
  desc <- var_desc[i]
  
  model <- lm_robust(as.formula(paste(v, "~ Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu")), data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster = cluster, weights = webal)
  
  models[[desc]] <- model
}

plot11 <- modelplot(models, coef_omit = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)") +
  labs(y = 'Treatment',
       title = "Attitudes toward the EU", subtitle = "EU Level") + geom_vline(xintercept = 0, color = 'grey')+
  coord_flip() + 
  theme(legend.position = "right", legend.justification="left") +
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        text = element_text(size = 12))+
  labs(x = NULL)
plot11
ggsave("out/other_outcomes.eps", plot = plot11, device = "eps")

########################################### Table 14 ###########################################
dat$agriculture <- NA
dat$agriculture[dat$qa19_1 == 1] <- 4
dat$agriculture[dat$qa19_1 == 2] <- 3
dat$agriculture[dat$qa19_1 == 3] <- 2
dat$agriculture[dat$qa19_1 == 4] <- 1
dat$agriculture <- as.numeric(dat$agriculture)

control_group <- dat[dat$Treatment == 0, ]
if (all(is.na(control_group))) {
  stop("Control group contains only NAs")
}

standardize <- function(var) {
  mean_var <- mean(control_group[[var]], na.rm = TRUE)
  sd_var <- sd(control_group[[var]], na.rm = TRUE)
  
  if (is.na(mean_var) || is.na(sd_var) || sd_var == 0) {
    warning(paste("Mean or SD is NA or SD is 0 for variable", var))
    return(rep(NA, nrow(dat)))
  }
  
  return((dat[[var]] - mean_var) / sd_var)
}

dat$agriculture_st <- standardize("agriculture")

modelpb_basic<-lm_robust(agriculture_st~Treatment, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_extended<-lm_robust(agriculture_st~Treatment + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_full<-lm_robust(agriculture_st~Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_balanced<-lm_robust(agriculture_st~Treatment + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster, weights = webal)


texreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Treatment and agriculture as cause of climate change", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-placebo",
       float.pos = "h",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

screenreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Treatment and agriculture as cause of climate change", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-placebo",
       float.pos = "h",
       file = "out/Table14.txt",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

########################################### Table 15 ###########################################
dat$Alternative_Treatment1=0
dat$Alternative_Treatment1[dat$Date < as.Date("2022-03-03")]=0
dat$Alternative_Treatment1[dat$Date == as.Date("2022-03-03")]=0
dat$Alternative_Treatment1[dat$Date > as.Date("2022-03-03")]=1

modelpb_basic<-lm_robust(eu_support_st~Alternative_Treatment1, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_extended<-lm_robust(eu_support_st~Alternative_Treatment1 + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_full<-lm_robust(eu_support_st~Alternative_Treatment1 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_balanced<-lm_robust(eu_support_st~Alternative_Treatment1 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster, weights = webal)


texreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Fake treatment and EU support", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-fake",
       float.pos = "h",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

screenreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Fake treatment and EU support", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-fake",
       float.pos = "h",
       file = "out/Table15.txt",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

########################################### Table 16 ###########################################

dat$Treatment3=0
dat$Treatment3[dat$Date > as.Date("2022-02-23")]=1
dat$Treatment3[dat$Date == as.Date("2022-02-24")]=NA

modelpb_basic<-lm_robust(eu_support_st~Treatment3, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_extended<-lm_robust(eu_support_st~Treatment3 + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_full<-lm_robust(eu_support_st~Treatment3 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_balanced<-lm_robust(eu_support_st~Treatment3 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster, weights = webal)

texreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Treatment effect when 24 February is removed and EU support", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-24rem",
       float.pos = "h",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

screenreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Treatment effect when 24 February is removed and EU support", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-24rem",
       float.pos = "h",
       file = "out/Table16.txt",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

########################################### Table 17 ###########################################

dat$Treatment2=0
dat$Treatment2[dat$Date > as.Date("2022-02-24")]=1

modelpb_basic<-lm_robust(eu_support_st~Treatment2, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_extended<-lm_robust(eu_support_st~Treatment2 + day, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_full<-lm_robust(eu_support_st~Treatment2 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster)
modelpb_balanced<-lm_robust(eu_support_st~Treatment2 + day + age + female + did_not_attend_school + primary_school + sec_grades_6_7 + sec_grades_1_5 + national_diploma + higher_national_diploma + bachelors_degree + masters_degree + doctoral_degree + education_level_10 + education_level_11 + education_level_12 + education_level_13 + education_level_14 + education_level_15 + education_level_16 + education_level_17 + education_level_18 + other_education + self_employed + managers + other_white_collar + manual_workers + house_persons + unemployed + retired + students + married_remarried + married_single_partner + single + divorced_separated + widow + other_marital_status + rural + small_town + large_town + internetuse + non_eu, data=dat, `se_type` = "stata", fixed_effects = ~ factor(country_code), cluster=cluster, weights = webal)

texreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Treatment on 25 February and EU support", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-25feb",
       float.pos = "h",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

screenreg(list(modelpb_basic, modelpb_extended, modelpb_full, modelpb_balanced), booktabs = TRUE, dcolumn = TRUE, include.ci = FALSE,
       caption = "Treatment on 25 February and EU support", stars = c(0.001, 0.01, 0.05), include.rmse=FALSE, 
       omit.coef = "day|age|female|did_not_attend_school|primary_school|sec_grades_6_7|sec_grades_1_5|national_diploma|higher_national_diploma|bachelors_degree|masters_degree|doctoral_degree|education_level_10|education_level_11|education_level_12|education_level_13|education_level_14|education_level_15|education_level_16|education_level_17|education_level_18|other_education|self_employed|managers|other_white_collar|manual_workers|house_persons|unemployed|retired|students|married_remarried|married_single_partner|single|divorced_separated|widow|other_marital_status|rural|small_town|large_town|internetuse|non_eu|(Intercept)", label = "tab:estimates-rc-25feb",
       float.pos = "h",
       file = "out/Table17.txt",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Time" = c("-", "Yes", "Yes", "Yes"),
                              "Covariates" = c("-", "-", "Yes", "Yes"),
                              "Weight" = c("-", "-", "-", "Yes")))

########################################### Table 18 ###########################################

bertelsmann_df = read_rds("~/Dropbox/Russia Ukraine EU/Replication/data/bertelsmann.rds")
bertelsmann_df$Remain=0
bertelsmann_df$Remain[bertelsmann_df$BM_EUREF=="2 - I would vote for my country to â€‹stay â€‹in the European Union"]=1
bertelsmann_df$EU_Active=0
bertelsmann_df$EU_Active[bertelsmann_df$BM_EUWORLD=="1 - I agree completely"]=1

reg0 = lm_robust(Remain ~ treatment, fixed_effects =country, data=bertelsmann_df, weights=survey_weight, se_type = "stata", subset = wave <= 26)
reg1 = lm_robust(Remain ~ treatment, fixed_effects =country, data=bertelsmann_df, weights=survey_weight, se_type = "stata", clusters = cluster, subset = wave <= 26)
reg2 = lm_robust(Remain ~ treatment+age + women + education + location, fixed_effects =country, weights= survey_weight, data=bertelsmann_df, subset = wave <= 26, clusters = cluster, se_type="stata")
reg3 = lm_robust(Remain ~ treatment+age + women + education + location, fixed_effects =country, weights = double_weights, data=bertelsmann_df, subset = wave <= 26, clusters = cluster, se_type="stata")

texreg(list(reg0, reg1, reg2, reg3), 
       include.ci = FALSE, digits = 3,
       caption = "Bertelsmann Dataset - Preferences to Remain in the EU",
       custom.model.names = c("Basic model", "Extended model", "Full model", "Balanced model"),
       omit.coef = c('date|age|women|education|location'),
       label = "table:bartel_remain",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Wave Cl. SE" = c("No", "Yes", "Yes", "Yes"),
                              "Covariates" = c("No", "No", "Yes", "Yes"),
                              "NN Weights" = c("No", "No", "No", "Yes")))

screenreg(list(reg0, reg1, reg2, reg3), 
       include.ci = FALSE, digits = 3,
       caption = "Bertelsmann Dataset - Preferences to Remain in the EU",
       custom.model.names = c("Basic model", "Extended model", "Full model", "Balanced model"),
       omit.coef = c('date|age|women|education|location'),
       label = "table:bartel_remain",
       file = "out/Table18.txt",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Wave Cl. SE" = c("No", "Yes", "Yes", "Yes"),
                              "Covariates" = c("No", "No", "Yes", "Yes"),
                              "NN Weights" = c("No", "No", "No", "Yes")))

########################################### Table 19 ###########################################

reg0 = lm_robust(EU_Active ~ treatment, fixed_effects =country, data=bertelsmann_df, weights=survey_weight, se_type = "stata", subset = wave <= 26)
reg1 = lm_robust(EU_Active ~ treatment, fixed_effects =country, data=bertelsmann_df, weights=survey_weight, se_type = "stata", clusters = cluster, subset = wave <= 26)
reg2 = lm_robust(EU_Active ~ treatment+age + women + education + location, fixed_effects =country, weights= survey_weight, data=bertelsmann_df, subset = wave <= 26, clusters = cluster, se_type="stata")
reg3 = lm_robust(EU_Active ~ treatment+age + women + education + location, fixed_effects =country, weights = double_weights, data=bertelsmann_df, subset = wave <= 26, clusters = cluster, se_type="stata")

texreg(list(reg0, reg1, reg2, reg3), 
       include.ci = FALSE,
       caption = "Bertelsmann Dataset - Should EU be more active in global affairs?",
       custom.model.names = c("Basic model", "Extended model", "Full model", "Balanced model"),
       omit.coef = c('date|age|women|education|location'),
       label = "table:bartel_remain",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Wave Cl. SE" = c("No", "Yes", "Yes", "Yes"),
                              "Covariates" = c("No", "No", "Yes", "Yes"),
                              "NN Weights" = c("No", "No", "No", "Yes")))

screenreg(list(reg0, reg1, reg2, reg3), 
       include.ci = FALSE,
       caption = "Bertelsmann Dataset - Should EU be more active in global affairs?",
       custom.model.names = c("Basic model", "Extended model", "Full model", "Balanced model"),
       omit.coef = c('date|age|women|education|location'),
       label = "table:bartel_remain",
       file = "out/Table19.txt",
       custom.gof.rows = list("Country FE" = c("Yes", "Yes", "Yes", "Yes"),
                              "Wave Cl. SE" = c("No", "Yes", "Yes", "Yes"),
                              "Covariates" = c("No", "No", "Yes", "Yes"),
                              "NN Weights" = c("No", "No", "No", "Yes")))
#

sink()
